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AN IMPROVED METHOD FOR SOLVING QUASILINEAR CONVECTION 
DIFFUSION PROBLEMS ON A COARSE MESH 

SARA POLLOCK 


Abstract. A method is developed for solving quasilinear convection diffusion prob¬ 
lems starting on a coarse mesh where the data and solution-dependent coefficients are 
unresolved, the problem is unstable and approximation properties do not hold. The 
Newton-like iterations of the solver are based on the framework of regularized pseudo¬ 
transient continuation where the proposed time integrator is a variation on the Newmark 
strategy, designed to introduce controllable numerical dissipation and to reduce the fluc¬ 
tuation between the iterates in the coarse mesh regime where the data is rough and the 
linearized problems are badly conditioned and possibly indefinite. An algorithm and up¬ 
dated marking strategy is presented to produce a stable sequence of iterates as boundary 
and internal layers in the data are captured by adaptive mesh partitioning. The method 
is suitable for use in an adaptive framework making use of local error indicators to de¬ 
termine mesh refinement and targeted regularization. Derivation and q-linear local con¬ 
vergence of the method is established, and numerical examples demonstrate the theory 
including the predicted rate of convergence of the iterations. 


1. Introduction 

This paper builds on the framework of ll2T]l and develops a nonlinear solver suitable 
for use in adaptive methods for quasilinear elliptic problems. The method is developed 
to stabilize the linearizations of nonlinear diffusion and convection diffusion problems, 
especially when there may be steep internal or boundary layers in the problem data. The 
sequence of linear problems encountered by a Newton-like method under these circum¬ 
stances takes the form of convection diffusion or reaction convection diffusion and the 
sequence of approximate solutions is subject to spikes, overshoots and spurious oscilla¬ 
tions in the convection-dominated regime. 

The present approach builds on a regularized version of the pseudo-transient contin¬ 
uation method as in for instance lUl [T9l IS [l2l |71 |2T]| and the references therein, and on 
each mesh refinement seeks a steady-state solution of the nonlinear evolution problem 
d/dt{Ru) + g{u{t)) = 0, for positive semidefinite regularization term R in the inter¬ 
est of solving g{u) = 0. In this analysis, further stabilization is introduced into the 
time discretization to address the problem of nonphysical oscillations. Time discretiza¬ 
tions featuring controllable high-frequency numerical dissipation are well known in the 
finite element analysis of structures as in for instance Il20l [T6l ITSl ITTI |4]], and a varia¬ 
tion related to these methods referred to here as the cr-split Newmark update is presently 
introduced. This method makes use of the controllable dissipation of the Newmark up¬ 
date and further controls fluctuations between the iterates by freezing a small fraction 
the linearization about a point with favorable properties. The cr-split Newmark method is 
derived, g-linear local convergence with a predictable rate is established, and the method 
is demonstrated on three variations of a model problem with steep internal layers. 
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The goal of the solver in the adaptive setting is to start on a coarse mesh where the 
problem data is generally not resolved and produce a sequence of transitional states 
which may not be accurate solutions to the coarse mesh problems, but do allow a poste¬ 
riori error indicators to detect the layers present in the problem data and refine the mesh 
to be globally fine enough for stability and locally fine enough to achieve accuracy and 
efficiency. The sequence of approximate problems can be classified into three phases, the 
initial phase where the mesh is globally too coarse and quadrature error is high and many 
features of the problem data remain undetected by the discretization, the pre-asymptotic 
phase where the mesh is fine enough for stability but the problem data is only partially re¬ 
solved, and the asymptotic regime where the standard existence, uniqueness and approx¬ 
imation properties hold. In the coarse mesh regime the solver uses as much stabilization 
as necessary to produce smooth transitional solutions; in the pre-asymptotic regime the 
solver adaptively reduces the added stabilization increasing both accuracy and the con¬ 
vergence rate as the data is resolved and the approximate problem becomes less rough; 
in the asymptotic phase the solver limits to a standard Newton method where the initial 
guess interpolated from the previous refinement is a good approximation to the solution 
and the iterations converge quadratically. 

The requirements of the adaptive method are that a local error indicator is computed 
on each refinement and a reasonable marking strategy in the sense of Il22ll is employed 
to determine the next mesh refinement. A modification of the standard adaptive marking 
strategy is proposed in which the marked set is determined by two parts: one that refines 
the elements with the largest indicators, and the other that refines a subset of the coarsest 
elements of the mesh when the residual from the final Newton-like iteration is significant. 
The method reported here improves on ll^ both in terms of efficiency and in terms of 
the strengths of near-singularities it is able to resolve. 

The remained of the paper is organized as follows. Section [3] shows the derivation of 
the pseudo-transient Newmark and cr-split Newmark methods, (12.81) and (12.91) . Section |4] 
demonstrates local g-linear convergence of the residual for these methods, both with rate 
q = 1 — 1 / 7 , with 7 the parameter from the Newmark update. Section |5] describes a 
basic algorithm to implement the solver in an adaptive method, and Section | 6 ] contains 
the results of numerical experiments using the described adaptive algorithm and (12.91) . 

The following notation is used in the remainder of the paper. The function g{u) refers 
to a specific problem or problem class and g{x) is used in the formal discussion of 
Newton-like methods. The nth iteration subordinate to the kih partition Tk is denoted 
x^, while is the nth iteration on a fixed partition and Xk is the final iteration on the 
kth mesh, taken as the approximate solution on Tk- In defining the weak and bilinear 
forms in the next section (u(x),v(x)) = u(x)v(x) dx, and similarly for vector-valued 
functions. 


2. Target problem class 

The nonlinear solver is developed to approximate solutions to the nonlinear problem 
g{u) = 0, for polygonal domain fl and g : X ^ Y* with g'{u) G C{X, Y*) for real Ba¬ 
nach spaces X and Y, particularly where g{u) takes the form of a quasilinear convection 
diffusion or diffusion problem in divergence form 

g{u) := —div(K(M)V«) + h{u) • Vu — f{x) = 0 in r 2 C m = 0 on dfl. ( 2 . 1 ) 
or 

g{u) := —div(K(M)V«) — f{x) = 0 in r 2 C m = 0 on dfl. ( 2 . 2 ) 
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Multiplication by a test function and integration by parts over the divergence term yields 
the weak form of ( 12 . Ih 

B{u, v) = {k{u)Vu, Vu) + {b{u) ■ Vu, v) for all V E ¥*, (2.3) 

and linearizing about uatwEX yields the bilinear form induced by g'{u) 

B'(u; w, v) = (K(M)Vtc, Vi;) + {k' (u)w'\/u, Vv) + {b{u) ■ Vic, v) + {b'{u)w'Vu, v). 

(2.4) 

For / G L 2 (r 2 ) nLoo(f2) and k{u) bounded away from zero with k(s), k,'{s) and k"{s) 
bounded as in |l3|, then (I2.2h has a unique solution u E with 2 < p < oo. 

Both (12. Ih and (12.21) fit into the context of Il24l with the assumption that k{u) is bounded 
and g'{u) : Hq{Q) is an isomorphism, in which case the solution u is an 

isolated solution. 

The discretized equation is, find Uh E such that B{uh, v) = 0 for all v E Vh where 
Xfi C X and 1/ C V are discrete finite element spaces with respect to triangulation 
where the family of triangulations {Th}o<h<i is regular and quasi-uniform in the sense 
of [|5]]. Existence, uniqueness and approximation properties of the discrete problems 
induced by (12.11) and (12.21) are found in Q, Il24ll and the references therein, assuming the 
mesh is fine enough. The problem is then to start on a coarse mesh; one that is not fine 
enough in terms of data resolution, stability or approximation properties, and build one 
that is. The goal of the solver is to navigate from a coarse to a sufficiently fine mesh 
where the approximation properties do hold, and to do so by building an both an efficient 
adaptive mesh and a reasonable initial guess to start the Newton-like iterations on each 
refinement along the way. 

The linearization of (12.21) has the form of a convection-diffusion equation, and (12.41) 
the linearized form of ( 12 . 11 ) has the form of a convection reaction diffusion equation 

a{w, v) + (3{w, v) + v) with 

a{w,v) \= {k{u)Xw,Vv) (2.5) 

/3(i(;, v) := {k\u)wXu, Vc) + {b(u) • Vic, v) ( 2 . 6 ) 

j{w,v) = {b'{u)wXu,v). (2.7) 

Using a standard Newton method to solve the nonlinear (12.3b with (12.4b generally does 
not work when there are steep layers present in the problem data, and coarse mesh ap¬ 
proximations of the problem are observed to be indefinite and ill-conditioned, consistent 
with the observations in [fT4]l . Many of the problems encountered including formation of 
spurious spikes, overshoots and instability are symptomatic of the corresponding linear 
convection-dominated problems |I3. Techniques form from the finite element analysis 
of structures ll20l[T6l[T5l[T7l l4ll use numerical integrators featuring high-frequency dissi¬ 
pation to capture lower frequency modes of the solution. The approach here investigates 
the use of the Newmark update and a stabilized variation of it as the time-integrators of 
a pseudo-transient continuation-like method as in lISlfTlIfil fT^I^ . 

The next section develops the two following methods to improve the convergence of 
the coarse-mesh iterates in the sequence of linearized problems. A positive semidefinite 
R is used to target specific degrees of freedom for regularization, and the role of R may 
be seen either from the homotopy perspective as a modification of the path between 
an initial and u* that solves g{u*) = 0, or from the regularization viewpoint as a 
penalty against certain characteristics of the iterates. In what follows, R is chosen based 
on the Laplacian to penalize against high curvature, with degrees of freedom selected 
for regularization based on an a posteriori error indicator. The selective approach to 
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regularization distinguishes the method presented here and in lllTII from the method of 
pseudo-transient continuation found elsewhere in the literature. For 7 > 1, a Newmark 
update generalizing a backward Euler discretization yields the iteration 

{anR + "f w'^ = —g{x^), = x'^ + w'^. ( 2 . 8 ) 

The proposed a-split Newmark update yields the iteration, 

{anR + "f {{1 — (j)g'{x) + ag'{x^)})w^ =—g{x'^), x^~^^ = x^ + w'^, (2.9) 

and analysis and demonstration of this last method is the focus of the remainder of the 
paper. Local convergence of (I2.8h is discussed along the way and local convergence 
of (12.9b is established by a perturbation of that result. The method based on (12.9b was 
also tested and found effective on a shifted p-Laplacian 

g{u) = —div { (c + IVm — (Vu — a)} — /, (2.10) 

and the related 

»(“)={('>+ ,+ivL«p ) ■ “>} ■ 

problems similar to those investigated in Q and ifT^ . which reside outside the current 
target class but still benefit from the regularization techniques described here when start¬ 
ing the adaptive method from a coarse mesh, especially if the coefficient e « 1 . 

3. Continuation methods 

The homotopy or pseudo-transient continuation method of stabilizing the Newton-like 
iterations for finding the solution x* of g{x) = 0 is developed by discretizing the ODE 

+ g{x{t)) = 0, x(0)=a:°, (3.1) 

at 

with a positive semidefinite linear functional R. In much of the literature, R is taken to 
be the identity or a scaled version thereof mini, and the references therein; however, the 
theory is developed for positive-definite functionals other than the identity Q where it is 
referred to as the “s” method, and positive semi-definite functionals in (^ET]]. This idea 
can be generalized to discretizing the ODE based on the normal equations formulation 
ofdO) 

d 

—R*Rx + g'{x{t))*g{x{t)) = 0, a:(0) = (3.2) 

with adjoint R* the formal adjoint of R and g'{x)* the formal adjoint of g'{x). As shown 
in II 2 T]], the discretization of (13.2b corresponds to the method of Tikhonov regularization. 

As discussed in for instance llll[8l[l9l|6]], letting x'^ a finite dimensional approxima¬ 
tion to x{tn) with Atn — tn+i — tn, the Standard method is to discretize (13.1b by a 
backward-Euler approximation to d/dt{R x) and a linearization of g'{x'^~^Y about x'^. In 
the case of discretizing (13.2b . p'(a:”+^)* is approximated by g'{x^)*. To increase stabil¬ 
ity of the linearized system, other discretizations of (13.2b are presently considered. The 
backward-Euler time discretization is replaced by the more general Newmark update, 
and the linearization of g{x'^~^Y is split about two distinct points, one the latest iterate 
and the other yielding a Jacobian with favorable properties. 

By the original method of backward-Euler discretization and first order Taylor expan¬ 
sion about the previous iterate, the resulting Newton-like iteration is given by 

(^-^R + g'ix^)^ = x^ + (3.3) 
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This method increases the stability of the linear system for positive definite and possibly 
semi-definite , but runs the risk of shifting the spectrum of the approximate Jacobian 
towards zero, making the condition significantly worse in the case that g'{u) is indefinite, 
resulting in large fluctuations between the iterates. To cope with this situation which 
occurs in the coarse mesh approximation of quasilinear problems, the iteration based on 
the normal equations enforcing the shift of the spectrum away from zero is given by 

( -^R*R + g'ix'^yg'ix'^U = -g'{x'^)g{x'^), x^+^ = x^ + w^. (3.4) 

This method is introduced in ifTil . where it is shown that (13.41) is also found by minimiz¬ 
ing the Tikhonov functional Ga{vo) for 

Ga{w) = \\g\x'^)w + an\\Rw\\l, (3.5) 

with an = and || • ||o the L 2 norm. As in IfTOll . the necessary and sufficient con¬ 

dition for the minimizer w is G'^{w) = 0, yielding (13.41) . While successfully increasing 
the stability of the system, this method remains unsatisfactory due to the increased com¬ 
plexity of solving the system based on the normal equations. A new method is now 
introduced which adds stability while preserving the sparsity of the system. The cost, 
as shown in Section |4] is trading the asymptotically quadratic convergence of (13.31) for 
g-linear convergence. The proposed algorithm of Section |5] updates the parameter 7 of 
the Newmark method on each mesh refinement as the sequence of linear systems stabi¬ 
lize until the method reduces to the original (13.31) for which quadratic convergence of the 
error is observed. 


3.1. Time discretization by the Newmark method. The Newmark method Il20ll dis¬ 
cretizes the time derivative ii = du/dthy 

-X^ = Mn {{I- 7 )t” + } . (3.6) 

For 7 = 1 , (13.61) reduces to the backward Euler discretization described above. As 
discussed in Il20ll this time integrator is second order accurate for 7 = 1 /2 and introduces 
nonphysical damping of the high frequency modes for 7 > 1/2 proportional to 7 — 
1 / 2 . More sophisticated a and 9 collocation methods as in |[T6l [l5l IHl incorporating the 
Newmark update are designed to further control numerical dissipation across targeted 
frequencies. In the current context the improved resolution and risk of overshoot of these 
methods designed with two time derivatives in mind is potentially of interest but their 
implementation in a pseudo-transient continuation setting is beyond the scope of this 
article. For stabilizing the transitional states of the sequence of coarse-mesh problems, 
damping of the high-frequency oscillations is the desirable property and more important 
than a higher order of accuracy. 

Solving (13.6b for 


x^+^ = 




^ — ^x^. 


7 

Applying Rx'^ = —g{xy and Rx^^^ = —g{x'^^y to (13.7b 

^ R{x^+^ - x^) = - g{x^^y 


'fAt 


7 


(3.7) 


(3.8) 


Linearizing g{x'^~^^) about x”, obtain the Newton-like iteration 


= x^ + 


(3.9) 
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For 7 > 1 this method comparatively scales down the influence of the right-hand side 
data; the danger is still allowing and arguably increasing the likeliness of shifting an 
indefinite left-hand side operator spectrum towards zero. By the same approach, the 
iteration based on the normal equations is given by 

=-g\x^)g{x^), . (3.10) 

Due to the overdamping effect of setting the parameter 7 > 1/2, the iterations (13.91) 
and (13.101) stabilize the solution in some situations where the iteration (I3.4h is infeasible. 
Similarly to (13.41) - (13.51) . the solution tc” of (13.101) is seen to be the minimizer of the 
generalized Tikhonov-type functional 

G'a(w) = \\g\x'^)w -g{x'^)\\l —\\Rw\\l, (3.11) 

7 7 

with an = 1 /Af„. 

3.2. The Newmark update with cr-splitting. The Newark discretization of the previ¬ 
ous section successfully introduces high-frequency dissipation increasing the stability of 
the linearized system but suffers the same drawback as the backward-Euler discretization 
requiring the formulation based on the normal equations in the case of a possibly indef¬ 
inite Jacobian to control highly unstable sequences of iterates. The cr-split Newmark 
update runs without the use of the normal equations: effectively freezing a small fraction 
of the Jacobian at a point with favorable properties dramatically reduces the fluctuations 
between iterates in the coarse mesh regime where boundary and internal layers are only 
partially resolved. 

As seen in the derivation, the method can be thought of as splitting the linearization of 
g'{x^^^) about two points, or more precisely as approximating Ax” by a combination of 
^n+i ^ ^j^g gj-gt jg ygg^j the backward-Euler update, the first two are 

used in the Newmark update with the second adding control of numerical dissipation; and 
the third is introduced here to reduce the fluctuation of the Jacobian outside the domain 
of convergence of the Newton-like iterations. The method is derived as follows, and is 
demonstrated in Section 0 to work without the use of the normal equations in situations 
where other methods including those involving the normal equations are seen to fail. 

Starting with the Newmark update (13.61) . each of the time derivative terms on the right 
is split into two parts 

- x” = At„{(l - 7 )(Tx” + (1 - 7)(1 - (j)x” + 7 ( 7 x”+^ + 7(1 - ct)x”+^} 

= At„{(l - 7 (j)x” + 7(1 - (j)(x”+^ - x”) + 7 ( 7 x”+^}. (3.12) 

Solving for and applying the relation = —g{x'^^^) 

—^RAx^ - - x”) = - ^(x”+^). (3.13) 

^aAtn cr 7(j 

Applying the relation Rx^ = —p(x”) on the left and linearizing g{x'^~^^) about x” on the 
right 

f^-i? + ^'(x”)') Ax” + ilA.^(p(x”+i) -p(x”)) = —p(x”). (3.14) 

V7crAt„ J a 70- 
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Linearizing both and g{x'^) about fixed x, g{x^~^^) — g{x‘^) = 

plying this and multiplying through by 'ja 


g'{x)Ax'^. Ap- 



(3.15) 


Equation (13.15b is the basis for the cr-split Newmark update. It contains four parame¬ 
ters: = 1 /Af„, 7 , cr and x', while some guidance is provided, a detailed analysis of 

these parameters will be addressed in subsequent work by the author. In the current re¬ 
sults x = 0 is used, and x = x^ the solution from the previous refinement interpolated 
onto the current mesh has also been observed to work; is chosen as in ll2T]l for the 
backward-Euler discretization, 7 > 1 may be set adaptively, increased to add stability 
then decreased to speed convergence; and cr > cto G (0,1) is adaptively sent towards one 
based on the norm of the latest residual on each Newton-like iteration. It is observed in 
numerical experiments that a should be close to one, and as seen in Section |4] this is a 
necessary for the asymptotic g-linear convergence at the rate g = 1 — I/ 7 , in agreement 
with the rate found using the Newmark discretization without the cr-splitting. 


4. Local convergence 


Local convergence of the residual is established for both algorithms (12.8b based on the 
Newmark update, and (12.9b based on the Nemark update with cr-splitting of the Jacobian. 
The second result follows as a perturbation of the first for cr close to unity, relying on 
a weaker set of assumptions. Both results require the same Lipschitz condition on the 
Jacobian. Denote the open ball about x by B{x,e) = {|/||| X — y\\ < e}. 

Assumption 4.1. There exist e > 0 50 that for all x,y G B{x*, e) 


\\ 9 \x) - g'{y)\\ < ujl\\x - y\\forallx,y G B{x\£). 


Assumption 4.2. (cf. Assumptions 2.2-2.3 o/|l6]], and Assumptions 4.1 and 4.4 of ll^ ). 
There exists {3 > b so that for positive semidefinite i?, 7 > 1, and for all 0 < < aM, 
then for all X G B{x*,£): 

1 ) anR + 'j g'(x) is invertible. 


2) \\(anR-\-'yg'ix)) ^\\ < M^. 

3) \\{anR){anR + 'yg'{x))-^\\ < 


— I^l3-f/an ‘ 


First, g-linear convergence of the residual is established with rate q = 1 — I/ 7 , for 
7 > 1. Linear and asymptotic quadratic convergence for the case 7 = 1 is shown 
in ll2T]l . following from convergence of the error. In the present discussion, convergence 
of the error is neither used nor shown: it is observed in numerical experiments that the 
residual decreases at the predicted rate over iterations where \\x^^^ — x'^\\ may not yet 
be decreasing; in contrast, for the case 7 = 1 the same quantity displays asymptoti¬ 
cally quadratic convergence to zero together with the residual. The following proof is a 
variation of Theorem 2.12 of |l7]|, where Assumptions 14. 1 1 and 14.21 replace the affine con- 
travariant Lipschitz condition used in that version of the Newton-Mysovskikh Theorem 
for the standard Newton method. 

Theorem 4.3. Let an < ||g(a:”)|| < a m and let Assumptions l4!7] and 14.21 hold. Define 
the open set S by 



(4.1) 
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and suppose G S; then by iteration (12.81) . G S, and the sequence of residuals 
converges q-linearly to zero with asymptotic rate q = 1 — l/j. 

Proof. Let Ax’* = x’*+^ — x**. Then iteration (12.81) is given by 

Ki? + 7 p'(x’*)) Ax’* =-p(x’*), x’*+^ = x” + Ax’*. (4.2) 

By the integral mean-value theorem 

g{x^ + AAx”) = g{x^) + [ g'{x^ + tAx**) Ax’* dt 


= £/(x’*) + A£/'(x’*) Ax’* + / (p'(x’* + tAx”) - £/'(x’*)) Ax’* dt 

Jo 

= fif(x’*) + a^R {ocnR + 

+ [ {g\x^ + tAx^)-g'{x^))Ax^dt (4.3) 

Jo 

Applying Assumption 14.21 (3) to the second term and the Lipschitz condition 14.11 to the 
third term of (14.3b 


|| 5 ((x’* AAx”)!! < 1 




A 


7(1 +^7/«n) 


\2 

ik^")ii + ^iiA^ir- 


(4.4) 


By the iteration (14.2b and Assumption l4Al 121. ||Ax|p < M^\\g{x‘^)\\‘^, yielding 


|| 5 ((x’* AAx**)!! < || 5 ((x’^ 


A 


\OLr^ 


7/ 7(«n + ^7) 


+ 






(4.5) 


By the assumption < || 5 f(x’*)||, and for A G [0,1] 

A 


||p(x’* + AAx’*)|| < \\g{x^‘ 


1 - 


7 




for each A G [0,1] such that x”+t Ax’* G S for all t G [0, A]. By the logic of (7]] Theorem 
2.12, proceed by contradiction and assume x’*+^ ^ S; then there is a smallest A G (0,1] 
with g{x‘^ + AAx**) G dS. For that A 

'A^ f 1 , 


||^(x’* + AAx’*)|| < ||^(x’*)|| 
< 






(4.7) 


a contradiction. This shows x’*"''^ G S. To establish the rate of convergence, set A = 1 
in l4.6[ 

1 \ yojLM^ ' 

,7/ V^7 ^ 2 


||^(x’*+i)|| < ||p(x’ 


1 - 


1 


7 


ll^(^”)ll , (4.8) 


which shows both that ||fi'(a:’*''‘^)|| < ||p(x’*)|| and the asymptotic g-linear rate of g = 
1 — 1/7 as ||g(a:’*)|| ^ 0 . □ 

Remark 4.4. It is observed that the method converges with the predicted rate when the 
iterates x^ f S as given by (14.1b . The lapse in the theory appears to be the bound on 
Ijx** II which apparently converges within a smaller set as compared to the residual. 
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The asymptotic convergence of the cr-split algorithm is addressed in the next theo¬ 
rem as a perturbation of Theorem 14.31 While it is not necessary in the proof for g'{x) 
to be positive definite, x is ideally chosen so that g'{x) improves the condition of the 
approximate Jacobian, allowing Assumptions 14.51 to hold with better constants than l4.2[ 

Assumption 4.5. (c.f. Assumption \4.2\) . There exist /3 > 0 i'c that for Q < uq < a < 1, 
fixed X, positive semidefinite R, 7 > and for all 0 < «„ < aM, then for all x e 

B{x*,e): 

1) UnR -|- 7{(1 — a)g'{x) — ag'{x)} is invertible. 

2) ||(Q;„i? + 7 {(l - a)g'{x) - ag'[x)})-^\\ < 

3) \\(anR)ianR + 7{1 - - (Tg'{x)})-^\\ < 

As with Assumption 14.21 (3), the third clause agrees with the similar stability bound 
of iU] Assumption 2.3, with replaced by an/T- Local convergence of the a-split 
Newmark algorithm is established with the same asymptotic rate as in the previous result. 


Theorem 4.6. Let an < || 5 f(a:”)|| < um and let Assumptions 14. 1 1 and 14.51 hold. Define a 
by 


a = max 



Ko J ’ 


(4.9) 


for a given Kq. Then there exists > 0 such that for x^ in the open set S given by 


S 


|x G B{x*,s) 


||c/(x)|| < , 


(4.10) 


and defined by iteration (12.91) . it hold that G S, and the sequence of residuals 
converges q-linearly to zero with asymptotic rate g = 1 — I/ 7 . 


Proof. Let — x'^. Then iteration (12.91) is given by 

(a^i? + 7{(1 — (j)g'(x) + (jg'(x”')}) Ax”'= —g(x’*), x”'*'^ = x” + Ax”. (4.11) 

Much of the proof parallels Theorem 14.31 and is summarized here. Starting with the 
integral mean-value theorem 


g(x” -|- AAx”) = (/(t”) -|- f g'{x‘^ + tAx”)Ax” dt 

Jo 


— g{x^) -l- Ag'(x”)Ax” -f / {g'{x'^ + fAx”) — g'(x”)) Ax” dt. 

Jo 

(4.12) 

By iteration (14.1 Ih 

g'(x”)Ax” = —^g(x”)-^a„i?Ax” — -— —g'{x)Ax‘^. (4.13) 

ya ya a 

Bounding the second term on the right of (14.131) by Assumptions |43] (3) 

\— 1 /_nA 


|a„i?Ax”|| = \\anR{anR + y{{l — (T)g\x) + (7 gfx'^)}) g(x^ 

1 


< 


1 + (3y! a, 




(4.14) 
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and by As sumptions 14.11 and 14.51 (2), then applying < || 5 '(x^)|| and A G [0,1] 

||,(." + AA.»)|| < m { (l - ^) + + (^) 


+ 




^CF 


M^nw 


< ||^(:r’^)|| { (l - ^) + 

+ (^:^) + —^^^ll^(^”)ll} • (4.15) 

Applying (14.91) to the quantity (1 — a)/a under the assumption 1 — \\g{x'^)\\/K q > ao 
and expanding in orders of ||p(x”)||/Ao < 1, 


(T Ko {l-\\g{x-)\\/Ko) Ko \ \ Ko )) 

M^nw , ^( Ux-)\\ V 

Ko \ Ko ) ' 

Similarly 

1 _ ^ ^ + (9 (^MGMV 

cr Kq \ Kq ) 

Applying (I4.16h and (I4.17h to (I4.15h . then for any fixed P > 3 


(4.16) 


(4.17) 


||^(a;" + AAx")||<||^(a;")|| 

'Koih^iXLMl/2 + 1 ) + ^-i{M,4g{x)\\ - 1 ) 


i--')+A||«K)II 

7/ Pi 


xP 


Kof^l 


PO{Ux-)\\^) 


(4.18) 


Supposing 

II ^ n^ii / r = Ko^jP 

Ko{f3^‘^u^Mll2p\)PM{M,^\g{x)\\-\Y 

it follows that 

||9(i" + AAi")|| < ||j(,T")|| (l - + 0(||j(x“)|p) , 

SO there exists (ii G (0, (io] for which ||fi'(a;”)|| < ()i implies 

+ AAi")|| < ||<,(i")|| (^1 - . 

The result follows by the same logic as in the proof of Theorem 14.31 


(4.19) 


(4.20) 


(4.21) 

□ 


Remark 4.7. It is observed in some problems that the residual || 5 f(x’*)|| can decrease at 
a rate slightly better than the one predicted when < 77 ^!, whereas for a = 1 the rate is 
as predicted. 
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5. Algorithm 

The solver using the cr-Newmark iteration (12.91) may be implemented in an adap¬ 
tive method according to the following basic algorithm. To exploit both the stability 
of method (12.91) and the quadratic convergence of (13.3h . the parameters 7 and R are up¬ 
dated on each refinement, and the parameters a and a are updated on each iteration of 
the solver. The main steps of the adaptive method are as follows with the exit criteria 
and parameter updates specified below. An example of a targeted regularization term R 
based on a posteriori error indicators is given in Section | 6 l 

Algorithm 5.1 (Basic algorithm for (I2.9M . Start with initial x^, 7 , (Tq. On partition 
Tk, A: = 0, 1 , 2 ,... 

1) Compute R, g'{x). 

2) SetaQ = |b(ic°)||. 

3) While Exit criteria \572\ are not met on iteration n — 1: 

(i) Set cr according to (14.91) . 

(ii) Solve {driR + + (1 — a)g'{x)})Ax"' = g{x^). 

(Hi) Update = x” + Ax”. 

(iv) Update an- 

4) Update 7/or partition Tk+i according to (15.3h 

Criteria 5.2 (Exit criteria). Given a user set tolerance tol, an accepted rate of conver¬ 
gence given by Qacc = 1 “ V {^l) for some constant M, e.g., M = 2, and a maximum 
number of iterations either chosen as a constant or based on the predicted rate of con¬ 
vergence, exit the solver on partition Tk after calculating iterate when one of the 

conditions holds. 

1) || 5 ((x”+^)|| < tol. 

2) (i) |b(x”+i)|| < \\g{x^)\\,AND 

(ii) |ifif(x”)|| < min{|| 5 f(x°)||, || 5 ((xfc_i)||}, AAD 

(Hi) \\g{x^+^')\\/\\g{x^)\\<qacc,AND 

(iv) ii^(x”+^)ij/jj^(x”)jj > ||^(x”)||/||^(x”-^)||. 

3) Maximum number of iterations exceeded. 

In terms of the three phases of the solution process in ll^ . the final asymptotic regime 
is characterized by the iterations terminating by Criteria l5^ (1): iterations in the pre- 
asympototic regime terminate with Criteria 15^ (2): and iterations terminate with a mix 
of Criteria |5]2] (2) and (3) in the initial phase. 

The second exit criterion 15121 (^21 merits some explanation, as it allows the iteration to 
end once reduction of the residual has slowed indicating the iterate has a attained a rea¬ 
sonably stable configuration from which a good prediction about where to refine the mesh 
may be determined. The first two clauses require the residual is decreasing, and has de¬ 
creased below the level given by the previous iterate with at least as much decrease from 
the initial iteration on the current partition and the residual from the previous partition, if 
it is well defined. The third criterion requires error reduction at or close to the predicted 
rate, and the last that the error reduction between iterates is slowing down. These four 
criteria together assure the sequence of transitional states is not getting further in the 
sense of the solver’s residual from a converged solution, and prevent situations such as 
spikes propagating indefinitely across a sequence of partitions. While spikes, overshoots 
or other undesirable characteristics may propagate through several refinements, such so¬ 
lutions will eventually not reduce the residual at the specified rate. When the adaptive 
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mesh is fine enough, these characteristics are observed to smooth out, and otherwise the 
solver eventually fails and the solution is reset and started again on a finer mesh. 

One of the main improvements of (12.91) over the regularized method of ll^ is the 
generation of more stable sequences of transitional iterates through the pre-asympototic 
phase. The sequence of pre-asymptotic approximate problems only partially resolve the 
data and as internal layers are uncovered over several refinements the condition of these 
problems is generally bad. The high-frequency dissipation of the Newmark strategy com¬ 
bined with the stabilization of the iterates produced by the cr-splitting allow sequences 
of solutions to propagate through this regime and the mesh to be marked via error indi¬ 
cators leading to the accurate solution of the problem on otherwise coarser meshes than 
possible if starting the solver on a mesh that is uniformly fine enough to resolve the data. 
In order to make use of the stability and dissipation properties of larger values of 7 = 7 ^ 
as well as the asymptotically quadratic convergence if the iterations are stable for 7 = 1 , 
the following minimal guidelines are presented, based on the termination criteria above. 

Update 5.3 (Newmark parameter 7 ). . Generally, if the predicted error reduction rate is 
achieved, 7 should be reduced; and if the iteration fails, more stability is needed and 7 
should be increased. 

(1) Exit criterium \5.2\ ( 1): Ify^ > 1 set 1 < ^k+i < Ik- 

(2) Exit criterium 15.21 12).• fjf ||5'(a:"^'''^^||/||fi'(a;”)|| is within tolerance ofl — 1 / 7 ^, set 
1 < Ik+i < Ik- 

(3) Exit criterium \5.2\ (3): Setyk^i > 7^. 

In practice, the residual tends to get reduced below tolerance once 7 = 1 . In the 
results of Section 7 is decreased by two when the target rate is achieved, decreased by 
one if a stable rate lower than predicted is achieved, increased by one when the solver 
fails, and by two if the maximum number of iteration is exceeded while the iterations are 
converging below the acceptable rate. An initial 70 should be chosen large enough to see 
error reduction on the initial mesh, and not significantly larger. 

The regularization parameter «„ = 1 / is updated by the method described in || 2 T]], 
repeated here for convenience. In accordance with the convergence Theorems l4.3l and l431 
this strategy assures < ||(y((a:”)|| so long as the residual is decreasing. 

Update 5.4 (Tikhonov regularization parameter a). Set Pq = 1. For n > 1, 

= Pn\\9{x'")\\, with Pn = • 

\\g{x^ 1)11 

To reduce rapid fluctuation of Pn, correct to ensure Pn-i/2 < /3„ < 1 in the case that 
\\g{x^)\\ < \\g{x^-'^)\\ andPn < 2/3„_i if\\g{x'^)\\ > \\g{x^-'^)\\. 

5.1. Marking strategy. An a posteriori error indicator pt, T e Tk is assumed avail¬ 
able to determine adaptive mesh refinement and as one option for determining a targeted 
regularization term R. In the current results, standard local residual-based element indi¬ 
cators as in for instance |l23l E]] are used for both of these purposes, further described 
in Section m Other approaches to solver- and problem-specific regularization and mesh 
refinement are currently under investigation by the author. 

The goal of the marking strategy is to build a mesh that is globally fine enough for 
stability, and locally as fine as necessary to achieve the desired accuracy. To improve the 
efficiency of the method by increasing the stability of the transitional states in the pre- 
asymptotic phase, the following marking strategy is presented. Based on exit criteria 15^ 
there are three possible outcomes of the nonlinear solve on refinement Tk- 
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1) Exit criterium l5.2l tli: The iterate has residual ||p(xfc)|| < tol. 

2) Exit criterium l5.2l (2): The iterate has residual ||p(xfc)|| > tol. 

3) Exit criterium 15.21 (31: The solver failed and Xk is reset to zero. 

In the case of 15.21 (31. the coarsest set of elements in the mesh are refined. Failure of 
the solver reflects a global rather than a local problem. Starting on a sufficiently coarse 
mesh, several resets are expected. 

In the case of 15.21 (21 the error indicators rjT are computed, and the mesh is refined 
both according to the elements with the largest local indicators, and according to the 
coarsest elements with the largest local indicators. This strategy allows local refinement 
to take place in order to capture boundary and internal layers to attain eventual accuracy 
and efficiency while also building the adaptive mesh fine enough to achieve stability. 
In meshes that are too coarse to resolve the data and where the approximate problem 
may have coefficients based on highly inaccurate approximate solutions, nonphysical 
overshoots often develop in the iterates; while the error indicators in the vicinity of these 
spikes may be high, refining primarily in these regions exacerbates the problem. As in 
case 15.21 (3'). a large residual || 5 f(ufc)|| predicts a global issue with the mesh. However, 
valuable information about the near-singularities in the problem data can be predicted 
from the non-converged iterates, so some local refinement can build a more efficient 
mesh. 

In the case of 15.21 (II the error indicators are computed, and the mesh is refined with re¬ 
spect to the largest local indicators. Any reasonable marking procedure Il22ll for cases l5.2l 
(l)-(2) may be applied; in particular for case 15.21 (21 both the element with the largest lo¬ 
cal indicator must be marked, as well as the coarsest element with the largest local indi¬ 
cator. In the current results, the Dorfler marking strategy is used with parameter 9, which 
for case 15.21 (21 is split into 9 = Oc + Op and the marked sets Aip Tk and Me Tk 
are chosen by sets of least cardinality with 


Vt>^fY Y 9t>^cY 

TeMp TeTk TeMc TeTk 

and the marked set Ad = Adi? U Adc- In case 15.21 (1'). 9f = 9 and 9c = 0. An heuristic 
choice of ^{9) = 9c is given in Section |6l 


6 . Numerical Examples 

The nonlinear solver (12.91) is demonstrated on three problems with different struc¬ 
ture in their internal layers. In Example 16.21 results are reported for a model nonlinear 
convection-diffusion problem of the form g{u) -0 which has a smooth sinusoidal solu¬ 
tion. The first variation on the model problem, Example lOl demonstrates the algorithm 
on the same differential operator with the problem data chosen to generate a higher fre¬ 
quency solution. Example l6^ shows the results for a related nonlinear diffusion problem 
with two concentric internal layers. The adaptive finite element method is implemented 
using the finite element library FETK IfTSlI and a direct solver is used on each linear sys¬ 
tem. Both trial and test spaces are taken as the linear finite element space 14 consisting 
of Lagrange finite elements Pi over partition Tk that satisfy the homogeneous Dirichlet 
boundary conditions. 
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The local a posteriori residual-based indicator, and corresponding jump-based indica¬ 
tor for element T E Tk with Ht the element diameter are given by 

Cr(^) = ^T||^r(^)||i2(9T) (6-1) 

Vri'^) — ^t\\9{'^)\\‘l2(t) + Ct(^)) (6-2) 

Jt{v) [[[k(i;)Vi; • nJgT, with jump [[0]aT := limj^o + tti) — 4>{x — tn), where n 
is the appropriate outward normal defined on dT. The error estimator on partition Tk is 
given by the I 2 sum of indicators = YhxeTk similarly for 

On each mesh partition Tk the penalty matrix R = Rk is a localized version of the 
Laplacian stiffness matrix denoted Let a diagonal matrix of zeros and ones, 

Vj a vertex of subordinate to partition Tk, and Cr = Criu^)- Then set 

4 = \/media„,.,r.(Cr), and A = | V • 
and 

nioc _ / 0) if Ct < i^k for each element that contains vj as a vertex, , 

R ^ 1, otherwise. ^ ‘ ^ 

Then set Rk = D^oc j^giobai p)ioc^ 


Remark 6.1. In lllTII a similar strategy is applied using the full indicator r}T as opposed 
to the C,T the jump-term alone. This modification was made to select degrees of freedom to 
penalize against curvature that show the greatest disparity in curvature as predicted by 
the jump term. The role of the selection function ip is to regularize based on spikes in the 
indicators leaving the stable dof unregularized to speed convergence. As the algorithm 
approaches the asymptotic phase the selection process ip becomes unimportant because 

«n < lb(a;”)|| ^ 0. 


In these results the Dorfler parameters 9c -\- Op = 0 with ^ = 0.6, for refinement k + 1 
is determined according to 

9c = ^{9) = ^ Q -f ;^arctan(||fif(«fc) 11/100 - 7r/2)^ . (6.4) 

This form of 4>(0) is chosen to ensure a significant fraction of the marked set is from the 
coarse mesh until the residual drops below a given threshold. In the following examples 
the parameter for setting a by formula (14.9b are a > oq = 0.9 and Kq = 2000. The 
tolerance for the nonlinear solver is tol = 10“^. 


Example 6.2 (Nonlinear convection diffusion). Consider the nonlinear convection dif¬ 
fusion equation onVt = [0,1] x [0,1], 

g{u) := —div(K(M) V«) + h{u) • Tu — f{x, y) = 0 E fl, u = 0 on dfl. (6.5) 
The nonlinear diffusion coefficient is given by 

k { s ) = k A - 77-7-777, with a = 0.5, and k = 1 . (6.6) 

((e+ (s-o)2) 

The nonlinear convection term is given by 

b{s) = {{s-a),{s-a)Y, 

and load f{x,y) chosen so the exact solution u(x,y) = sin ttx sin iry. 


(6.7) 
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Existence of isolated solutions of (I6.5I) - (I6.7I) is discussed in Section |2] following Il24ll . 
assuming the mesh is fine enough. This problem without the convection term is investi¬ 
gated in ll2T]l . with e = 10 “^; the method here has been observed to run to convergence 
from a mesh of 144 elements with e = 8 x 10“^, and was tested on smaller values of e; 
however the size of the resulting meshes were unwieldy for test-problems. Data on the 
iterations in the pre-asymptotic regime is presented here for £ = 2 x 10 “^ starting on a 
mesh of 144 elements, with 7 = 10 and run until convergence; and for c = 6 x 10“^ run 
well into the asymptotic regime. In the pre-asymptotic case the iterations are terminated 
by exit criteria 151^ 121. and in most cases where 7 is reduced it is reduced by one as the 
iterates are converging at a rate better than predicted. On refinements 20,21 and 28 the 
final iterates are converging at or close to the predicted rate q = 1 — I/ 7 , and 7 is re¬ 
duced by two. Table [U summarizes the computation showing the final residual || 5 '(Mfc)||, 
the final ratio \\g{u^~^^) ||> the final value of a and the Newmark parameter 7 . 

Of note is the range in magnitude over which the residual shows monotonic and in 
fact predictable and stable decrease. On level 7, the first partition in the pre-asymptotic 
regime, the residual is reduced from an initial ||fi'(u°)|| = 1564.5, to a final|| 5 ((u®)|| = 
871.1, using a = (Tq = 0.9 on each iteration. The final pre-asymptotic solve reduces the 
residual from || 5 '(u°)|| = 4.6, to || 5 '(m^)|| = 1.4, with a > 0.996 and approaching one. 
On level 7 for which 7 = 20, the ratio of consecutive residuals stabilizes well below the 
predicted rate whereas in level 28, the predicted rate 1 — I/ 7 , with 7 = 3 , is achieved. 
Most of the pre-asymptotic solves exit after 3-4 iterations, and two of the solves take up 
to 9; generally the higher number of iterations corresponds to qualitative changes in the 
iterates such as the smoothing of spikes and changes in curvature illustrated for instance 
in the solution snapshots in Figure HI 

Setting £ = 6 X 10“^, a milder version of the same problem illustrates the reduction 
in the error from the initial through the pre-asymptotic and into the asymptotic regime 
which in this case runs from levels 2 to 19. The computation is again started on a mesh 
of 144 elements with 7 = 10. Figure [H shows a logarithmic plot of the error against 
the number of elements n compared with The increase in the error and estimator 

for the first several refinements demonstrates the sequences of partitions capturing the 
internal layer in the problem data; the first mesh in Figure [Hand corresponding solution 
in Figure |2j show the representative behavior as this occurs; the flatness at the top of 
the solution is characteristic of the penalization by the Laplacian against curvature. The 
second mesh of Figure [U and corresponding solution of Figure |3] show the approximate 
solution as the internal layer is better resolved, but spurious nonphysical oscillations are 
still apparent in the iterates; this mesh further illustrates the characteristic of the layer 
getting captured by the mesh, but not uniformly. The final mesh of Figure |2] and solution 
of Figure |3] are representative of the end of the asymptotic regime where the mesh now 
refines to capture the layer with increasingly uniform resolution and the solution has 
the overall correct shape. Within several iteration the approximate solution looks like a 
sinusoid and the Newton-like iterations solve to tolerance. It is noted in Figure (H that 
the error is stable for the first several iterations in the asymptotic regime while the error 
estimator reduces at a steady rate; presumably this occurs due to the presence of another 
nearby solution of (I6.5h . A similar phenomenon is observed in IIT^ in examples where 
the monotonicity is violated and the solution is known only to be locally unique. 

Example 6.3 (Higher frequency solution). A variation of Example 16.21 is considered 
where the problem is given by (I6.5I) - (I6.7I) . with e = 6 x 10“^, and the load f{x,y) 
chosen so the exact solution is u = sin(27ra:) sin(27r|/). 
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Level 

\\9M\\ 

||9(«"+l)|| 

O'k 

7fc 

7 

871.1 

0.90 

0.9 

20 

8 

786.8 

0.93 


20 

9 

647.0 

0.94 


19 

10 

535.0 

0.94 


18 

11 

445.0 

0.93 


17 

12 

368.3 

0.93 


16 

13 

305.2 

0.93 


15 

14 

275.7 

0.92 


14 

15 

231.3 

0.92 

0.9 

13 

16 

165.6 

0.91 

0.909 

12 

17 

137.7 

0.90 

0.924 

11 

18 

122.1 

0.90 

0.932 

10 

19 

98.4 

0.88 

0.944 

9 

20 

75.0 

0.87 

0.957 

8 

21 

59.0 

0.83 

0.964 

6 

22 

44.0 

0.75 

0.971 

4 

23 

23.5 

0.67 

0.983 

3 

24 

14.1 

0.90 

0.990 

3 

25 

7.7 

0.71 

0.995 

3 

26 

4.8 

0.68 

0.997 

3 

27 

2.5 

0.67 

0.998 

3 

28 

1.4 

0.67 

0.999 

3 

29 

1.9e-08 

8.1e-03 

1 

1 


Table 1. Summary of data from the pre-asymptotic regime for Exam- 
ple 16.21 with g = 2 x 10“^. 



Figure 1. error (below) and error estimator (above) against n 
where n is the number of elements, for Example 16.21 with £ = 6 x 10“^. 
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Figure 2. Mesh for Example 16.21 with e = 6 x 10 ^ after 6, 12 and 18 
adaptive refinements from an initial mesh with 144 elements. 



Figure 3. Snapshots of the finite element solution for Example |62] with 
£ = 6 X 10“^ in the pre-asymptotic regime after 6, 12 and 18 adaptive 
refinements from an initial mesh with 144 elements. 


Local uniqueness of the solution follows from ll^ . assuming the mesh is fine enough. 
This example features two disjoint internal layers in the problem data, one about each 
positive peak of the sinusoid. The decrease in the residual in the pre-asymptotic regime 
closely resembles the data presented in Tabled for Example l6.2[ For £ = 6 x 10“^, 
the pre-asymptotic phase for this example runs from levels 4 to 26 as opposed to 2 to 
19 for the sinusoid with a single peak; and the mesh in Example 16.31 maintains a maxi¬ 
mum meshsize of h = 8.6 x 10“^, whereas in Example 16.21 the asymptotic meshsize is 
h = 3.5 X 10“^, indicating two more coarse mesh refinements are taken to stabilize the 
iterates. This is to be expected due to the decrease in width of the corresponding lay¬ 
ers. Mesh partitions and snapshots of the pre-asymptotic iterates are shown in Figures |4] 
and S The characteristic flatness of the peaks due to the penalization against curvature 
is again observed in the first two iterates, shown respectively at the 8th and 12th adaptive 
refinements, and resolved by the 22nd, as seen on the right. The iterates show a qual¬ 
itatively different behavior than the nonphysical oscillations of Example 16.21 here, the 
downward peaks which should have the same magnitude as the upward peaks start as 
shallow and extend to their full depth as the refinements progress. 

Example 6.4 (Concentric layers). Consider the nonlinear diffusion problem on fl = 
[0,1] X [0,1] with two concentric internal layers. 


g{u) := —div(K(M) V«) — f{x, y) = 0 E fl, u = 0 on dfl, 


( 6 . 8 ) 
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Figure 4. Mesh for Example 16.31 with c = 6 x 10 ^ after 8, 14 and 22 
adaptive refinements from an initial mesh with 144 elements. 





Figure 5. Snapshots of the finite element solution for Example l6.3l with 
£ = 6 X 10“^ in the pre-asymptotic regime after 8, 14 and 22 adaptive 
refinements from an initial mesh with 144 elements. 


with the nonlinear dijfusion coefficient given by 
1 1 


k{s) = k + 


+ 


((6 +(.-a)2) ((e + (.-c)2) 


, with a 


The load f{x, y) is chosen so the exact solution u{x, y) = 


= 0.5, c = 0.8 and k 
sin(7rx) sin(7r|/). 


1 . 

(6.9) 


Local uniqueness of the solution again follows from ll24ll . assuming the mesh is fine 
enough. Similarly to the error decrease in the first example. Figure 0 shows a temporary 
leveling off of the error while the estimator decreases at a steady rate at the beginning of 
the asymptotic phase indicating the presence of a nearby solution. The discontinuity in 
the error estimator and corresponding jump in the error is due to the final coarse mesh 
refinement where the solution is reset at adaptive level 22, after which the mesh maintains 
a maximum meshsize of h = 1.7 x 10“^, half the meshsize of Example 16.21 Once in 
the pre-asymptotic regime the internal layers are progressively resolved and the residual 
drops below tolerance converging at a quadratic rate when 7 = 1, at adaptive level 37. 
Due to the aforementioned flattening of the iterates induced by the solver’s stabilization, 
the approximate solutions gradually increase in height as the regularization parameters 
are reduced; as such, the inner internal layer centered at u = 0.8, is uncovered later 
than the outer layer centered atu = 0.5. This phenomenon is illustrated in the adaptive 
meshes at levels 25, 30 and 35 shown in Figure |71 The corresponding iterates shown 
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Figure 6 . error (below) and error estimator (above) against n 
where n is the number of elements, for Example 16.41 with £ = 6 x 10“^. 



Figure 7. Mesh for Example lOI with e = 6 x 10 ^ after 25, 30 and 35 
adaptive refinements from an initial mesh with 144 elements. 





Figure 8. Snapshots of the finite element solution for Example 16.41 with 
£ = 6 X 10“^ in the pre-asymptotic regime after 25, 30 and 35 adaptive 
refinements from an initial mesh with 144 elements. 


in Figure [8] display the characteristic spikes between the outer layer and the boundary 
similar to those in Example 16.21 and in this case the peak of the sinusoid progressively 
resolves the curvature in two flatter regions. 
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7. Conclusion 

The (7-split Newmark method of regularized pseudo-transient continuation is intro¬ 
duced to solve the nonlinear convection diffusion problem g{u) = Q starting from a 
coarse mesh where the problem data is not yet resolved. The method is derived and local 
g-linear convergence in agreement with the observed rate is established. An algorithm is 
presented to fit the solver into a standard adaptive method, and is demonstrated on three 
variations of a problem with steep internal layers. The method including the adaptive 
update of the solver’s parameters on each refinement builds on the framework presented 
in [|2T]| . and is designed to effectively stabilize linearizations over rough problem data 
where spikes, overshoots and spurious oscillations otherwise prevent the sequence of 
transitional solutions from approaching an accurate approximation. The results here are 
an improvement first in terms of efficiency of the new solver which functions without 
the use of the normal equations; and second by an improved set of exit criteria for the 
solver which works together with an updated marking strategy to improve the stability 
of the method through the pre-asymptotic coarse mesh regime. Ultimately this yields a 
more efficient route to the asymptotic regime where the problem data is resolved and ap¬ 
proximation properties hold. Future work will address targeted local regularization and 
parameter selection of the nonlinear solver, and more general types of nonlinearites. 
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